{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Velocity Data Download Script for Validation\n",
    "#### Working Script to pull large velocity datasets and later compare to ICESAT2-derived velocities\n",
    "\n",
    "ICESat-2 hackweek  \n",
    "June 15, 2020  \n",
    "Lynn Kaluzienski"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import necessary modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os,re,h5py\n",
    "import requests\n",
    "import zipfile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import pyproj\n",
    "import scipy, sys, os, pyproj, glob\n",
    "import matplotlib.pyplot as plt\n",
    "from shapely.geometry import Point, Polygon\n",
    "# run matplotlib in 'widget' mode\n",
    "%matplotlib widget\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#! cd ..; [ -d pointCollection ] || git clone https://www.github.com/smithB/pointCollection.git\n",
    "sys.path.append(os.path.join(os.getcwd(), '..'))\n",
    "\n",
    "import pointCollection as pc\n",
    "#stuck here, not sure how to import point collection \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Magic function to enable interactive plotting (zoom/pan) in Jupyter notebook\n",
    "#If running locally, this would be `%matplotlib notebook`, but since we're using Juptyerlab, we use widget\n",
    "%matplotlib widget\n",
    "#%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "##Section where I played around with loading in .mat file, quickly abandoned\n",
    "#load in data from .mat file\n",
    "#Data originally in netcdf format, created subset in matlab for FIS region and saved as .mat file. \n",
    "import scipy.io as spio\n",
    "data_root='/srv/surface_velocity/contributors/lynn_kaluzienski/Data_test'\n",
    "#('Measures2.mat',squeeze_me=True)\n",
    "Measures2_Files= spio.loadmat('./Data_test/Measures2.mat')\n",
    "x=Measures2_Files['posx']\n",
    "y=Measures2_Files['posy']\n",
    "vx=Measures2_Files['vx_FIS']\n",
    "vy=Measures2_Files['vy_FIS']\n",
    "vel=Measures2_Files['vel_FIS']\n",
    "\n",
    "#xy=np.array(pyproj.Proj(3031)(x, y))\n",
    "plt.figure()\n",
    "ax.imshow(vel, extent=[0, 1, 0, 1])\n",
    "''';"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-86 -81 -81 -86 -86]\n",
      "[-55 -55 -65 -65 -55]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7ccce4c822e7407f8e08d7e76b17a95e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'cmap': 'viridis', 'clim': [0, 500], 'extent': array([-887878.65809562, -400567.81840096,  200333.90920048,\n",
      "        561654.87344315]), 'origin': 'lower'}\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Measures2 Velocity Map')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#From Ben Smith's code loading in .tif file, running into issues likely with directories\n",
    "data_root='/srv/shared/surface_velocity/FIS_Velocity/'\n",
    "#spatial_extent = np.array([-102, -76, -98, -74.5])\n",
    "spatial_extent = np.array([-65, -86, -55, -81])\n",
    "lat=spatial_extent[[1, 3, 3, 1, 1]]\n",
    "lon=spatial_extent[[2, 2, 0, 0, 2]]\n",
    "print(lat)\n",
    "print(lon)\n",
    "# project the coordinates to Antarctic polar stereographic\n",
    "xy=np.array(pyproj.Proj(3031)(lon, lat))\n",
    "# get the bounds of the projected coordinates \n",
    "XR=[np.nanmin(xy[0,:]), np.nanmax(xy[0,:])]\n",
    "YR=[np.nanmin(xy[1,:]), np.nanmax(xy[1,:])]\n",
    "#Originally tried to load data from a local directory, should change to shared directory\n",
    "Measures_vel=pc.grid.data().from_geotif(os.path.join(data_root,'Measures2_FIS_Vel.tif'), bounds=[XR, YR])\n",
    "\n",
    "# show the velocity map:\n",
    "plt.figure()\n",
    "Measures_vel.show(cmap='viridis', clim=[0,500])\n",
    "plt.plot(xy[0,:], xy[1,:])\n",
    "plt.title('Measures2 Velocity Map')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def atl06_to_dict(filename, beam, field_dict=None, index=None, epsg=None):\n",
    "    \"\"\"\n",
    "        Read selected datasets from an ATL06 file\n",
    "\n",
    "        Input arguments:\n",
    "            filename: ATl06 file to read\n",
    "            beam: a string specifying which beam is to be read (ex: gt1l, gt1r, gt2l, etc)\n",
    "            field_dict: A dictinary describing the fields to be read\n",
    "                    keys give the group names to be read, \n",
    "                    entries are lists of datasets within the groups\n",
    "            index: which entries in each field to read\n",
    "            epsg: an EPSG code specifying a projection (see www.epsg.org).  Good choices are:\n",
    "                for Greenland, 3413 (polar stereographic projection, with Greenland along the Y axis)\n",
    "                for Antarctica, 3031 (polar stereographic projection, centered on the Pouth Pole)\n",
    "        Output argument:\n",
    "            D6: dictionary containing ATL06 data.  Each dataset in \n",
    "                dataset_dict has its own entry in D6.  Each dataset \n",
    "                in D6 contains a numpy array containing the \n",
    "                data\n",
    "    \"\"\"\n",
    "    if field_dict is None:\n",
    "        field_dict={None:['latitude','longitude','h_li', 'atl06_quality_summary'],\\\n",
    "                    'ground_track':['x_atc','y_atc'],\\\n",
    "                    'fit_statistics':['dh_fit_dx', 'dh_fit_dy']}\n",
    "    D={}\n",
    "    file_re=re.compile('ATL06_(?P<date>\\d+)_(?P<rgt>\\d\\d\\d\\d)(?P<cycle>\\d\\d)(?P<region>\\d\\d)_(?P<release>\\d\\d\\d)_(?P<version>\\d\\d).h5')\n",
    "    with h5py.File(filename,'r') as h5f:\n",
    "        for key in field_dict:\n",
    "            for ds in field_dict[key]:\n",
    "                if key is not None:\n",
    "                    ds_name=beam+'/land_ice_segments/'+key+'/'+ds\n",
    "                else:\n",
    "                    ds_name=beam+'/land_ice_segments/'+ds\n",
    "                if index is not None:\n",
    "                    D[ds]=np.array(h5f[ds_name][index])\n",
    "                else:\n",
    "                    D[ds]=np.array(h5f[ds_name])\n",
    "                if '_FillValue' in h5f[ds_name].attrs:\n",
    "                    bad_vals=D[ds]==h5f[ds_name].attrs['_FillValue']\n",
    "                    D[ds]=D[ds].astype(float)\n",
    "                    D[ds][bad_vals]=np.NaN\n",
    "    if epsg is not None:\n",
    "        xy=np.array(pyproj.proj.Proj(epsg)(D['longitude'], D['latitude']))\n",
    "        D['x']=xy[0,:].reshape(D['latitude'].shape)\n",
    "        D['y']=xy[1,:].reshape(D['latitude'].shape)\n",
    "    temp=file_re.search(filename)\n",
    "    D['rgt']=int(temp['rgt'])\n",
    "    D['cycle']=int(temp['cycle'])\n",
    "    D['beam']=beam\n",
    "    return D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0744244a241a438285ded062dbd9dfa3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'cmap': 'viridis', 'clim': [0, 400], 'extent': array([-887878.65809562, -400567.81840096,  200333.90920048,\n",
      "        561654.87344315]), 'origin': 'lower'}\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fcc455ebcd0>]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load a line and plot\n",
    "\n",
    "data_root='/srv/shared/surface_velocity/'\n",
    "field_dict={None:['delta_time','latitude','longitude','h_li', 'atl06_quality_summary'],\\\n",
    "                    'ground_track':['x_atc','y_atc'],\\\n",
    "                    'fit_statistics':['dh_fit_dx', 'dh_fit_dy']}\n",
    "\n",
    "rgt = \"0848\"\n",
    "cycle=\"03\"\n",
    "filename = glob.glob(os.path.join(data_root, 'FIS_ATL06', f'*ATL06_*_{rgt}{cycle}*_003*.h5'))\n",
    "\n",
    "D=atl06_to_dict(filename[0],'/gt2l', field_dict=field_dict, index=None, epsg=3031)\n",
    "\n",
    "# show the velocity map:\n",
    "plt.figure()\n",
    "Measures_vel.show(cmap='viridis', clim=[0,400])\n",
    "plt.plot(xy[0,:], xy[1,:],'k')\n",
    "plt.title('Measures2 Velocity Map')\n",
    "\n",
    "plt.plot(D['x'],D['y'],'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b3bc18f5b8f540bab195b84146f269c7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fcc2a583150>]"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Interpolate the Measures velocities along the line and plot (note that these are speeds for now)\n",
    "\n",
    "V = Measures_vel.interp(D['x'],D['y'], gridded=False)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(D['x_atc'],V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# extract data velocity points\n",
    "\n",
    "idea from\n",
    "https://www.earthdatascience.org/courses/use-data-open-source-python/spatial-data-applications/lidar-remote-sensing-uncertainty/extract-data-from-raster/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import rasterio as rio\n",
    "import geopandas as gpd\n",
    "\n",
    "\n",
    "from geopandas import GeoDataFrame\n",
    "from shapely.geometry import Point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Dictionary to Pandas df\n",
    "df = pd.DataFrame.from_dict(D)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/pyproj/crs/crs.py:53: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6\n",
      "  return _prepare_from_string(\" \".join(pjargs))\n"
     ]
    }
   ],
   "source": [
    "#Convert points to shape file or geopandas\n",
    "\n",
    "geometry = [Point(xy) for xy in zip(df.x, df.y)]\n",
    "df = df.drop(['x', 'y'], axis=1)\n",
    "crs = {'init': 'epsg:3031'}\n",
    "gdf = GeoDataFrame(df, crs=crs, geometry=geometry)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "#create shapefile\n",
    "\n",
    "gdf.to_file('test_line_is2.shp', driver='ESRI Shapefile')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>delta_time</th>\n",
       "      <th>latitude</th>\n",
       "      <th>longitude</th>\n",
       "      <th>h_li</th>\n",
       "      <th>atl06_quality_summary</th>\n",
       "      <th>x_atc</th>\n",
       "      <th>y_atc</th>\n",
       "      <th>dh_fit_dx</th>\n",
       "      <th>dh_fit_dy</th>\n",
       "      <th>rgt</th>\n",
       "      <th>cycle</th>\n",
       "      <th>beam</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>4.387628e+07</td>\n",
       "      <td>-81.000132</td>\n",
       "      <td>-55.500558</td>\n",
       "      <td>99.545357</td>\n",
       "      <td>0</td>\n",
       "      <td>2.912429e+07</td>\n",
       "      <td>36.549355</td>\n",
       "      <td>0.000697</td>\n",
       "      <td>-0.000276</td>\n",
       "      <td>848</td>\n",
       "      <td>3</td>\n",
       "      <td>/gt2l</td>\n",
       "      <td>POLYGON ((-807437.572 554952.385, -807437.764 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.387628e+07</td>\n",
       "      <td>-81.000305</td>\n",
       "      <td>-55.500826</td>\n",
       "      <td>99.541039</td>\n",
       "      <td>0</td>\n",
       "      <td>2.912431e+07</td>\n",
       "      <td>36.477020</td>\n",
       "      <td>-0.004370</td>\n",
       "      <td>-0.000480</td>\n",
       "      <td>848</td>\n",
       "      <td>3</td>\n",
       "      <td>/gt2l</td>\n",
       "      <td>POLYGON ((-807424.574 554937.880, -807424.766 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.387628e+07</td>\n",
       "      <td>-81.000478</td>\n",
       "      <td>-55.501095</td>\n",
       "      <td>99.530708</td>\n",
       "      <td>0</td>\n",
       "      <td>2.912433e+07</td>\n",
       "      <td>36.402977</td>\n",
       "      <td>0.002140</td>\n",
       "      <td>-0.000554</td>\n",
       "      <td>848</td>\n",
       "      <td>3</td>\n",
       "      <td>/gt2l</td>\n",
       "      <td>POLYGON ((-807411.571 554923.378, -807411.764 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.387628e+07</td>\n",
       "      <td>-81.000652</td>\n",
       "      <td>-55.501362</td>\n",
       "      <td>99.531052</td>\n",
       "      <td>0</td>\n",
       "      <td>2.912435e+07</td>\n",
       "      <td>36.344711</td>\n",
       "      <td>0.000911</td>\n",
       "      <td>-0.000570</td>\n",
       "      <td>848</td>\n",
       "      <td>3</td>\n",
       "      <td>/gt2l</td>\n",
       "      <td>POLYGON ((-807398.559 554908.886, -807398.752 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4.387628e+07</td>\n",
       "      <td>-81.000825</td>\n",
       "      <td>-55.501629</td>\n",
       "      <td>99.546165</td>\n",
       "      <td>0</td>\n",
       "      <td>2.912437e+07</td>\n",
       "      <td>36.305172</td>\n",
       "      <td>0.004251</td>\n",
       "      <td>-0.000422</td>\n",
       "      <td>848</td>\n",
       "      <td>3</td>\n",
       "      <td>/gt2l</td>\n",
       "      <td>POLYGON ((-807385.533 554894.405, -807385.726 ...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     delta_time   latitude  longitude       h_li  atl06_quality_summary  \\\n",
       "0  4.387628e+07 -81.000132 -55.500558  99.545357                      0   \n",
       "1  4.387628e+07 -81.000305 -55.500826  99.541039                      0   \n",
       "2  4.387628e+07 -81.000478 -55.501095  99.530708                      0   \n",
       "3  4.387628e+07 -81.000652 -55.501362  99.531052                      0   \n",
       "4  4.387628e+07 -81.000825 -55.501629  99.546165                      0   \n",
       "\n",
       "          x_atc      y_atc  dh_fit_dx  dh_fit_dy  rgt  cycle   beam  \\\n",
       "0  2.912429e+07  36.549355   0.000697  -0.000276  848      3  /gt2l   \n",
       "1  2.912431e+07  36.477020  -0.004370  -0.000480  848      3  /gt2l   \n",
       "2  2.912433e+07  36.402977   0.002140  -0.000554  848      3  /gt2l   \n",
       "3  2.912435e+07  36.344711   0.000911  -0.000570  848      3  /gt2l   \n",
       "4  2.912437e+07  36.305172   0.004251  -0.000422  848      3  /gt2l   \n",
       "\n",
       "                                            geometry  \n",
       "0  POLYGON ((-807437.572 554952.385, -807437.764 ...  \n",
       "1  POLYGON ((-807424.574 554937.880, -807424.766 ...  \n",
       "2  POLYGON ((-807411.571 554923.378, -807411.764 ...  \n",
       "3  POLYGON ((-807398.559 554908.886, -807398.752 ...  \n",
       "4  POLYGON ((-807385.533 554894.405, -807385.726 ...  "
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# bufffer over points\n",
    "\n",
    "gdf.buffer(20).to_file('test_line_is2.shp')\n",
    "\n",
    "is2_poly =gdf.copy()\n",
    "is2_poly[\"geometry\"] = gdf.geometry.buffer(40) # 40 meters buffer for zonal statitics\n",
    "is2_poly.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If the dir does not exist, create it\n",
    "output_path = os.path.join(\"poly\")\n",
    "\n",
    "if not os.path.isdir(output_path):\n",
    "    os.mkdir(output_path)\n",
    "\n",
    "# Export the buffered point layer as a shapefile to use in zonal stats\n",
    "plot_buffer_path = os.path.join(output_path, \"plot_buffer.shp\")\n",
    "is2_poly.to_file(plot_buffer_path)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rasterstats contains the zonalstatistics function \n",
    "# that you will use to extract raster values\n",
    "import rasterstats as rs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract zonal stats\n",
    "Measures_vel_2=rio.open(os.path.join(data_root,'Measures2_FIS_Vel.tif'))\n",
    "\n",
    "MeASUREs_vel2 = rs.zonal_stats(plot_buffer_path,\n",
    "                                   SJER_chm_data,\n",
    "                                   nodata=-999,\n",
    "                                   affine=sjer_chm_meta['transform'],\n",
    "                                   geojson_out=True,\n",
    "                                   copy_properties=True,\n",
    "                                   stats=\"count min mean max median\")\n",
    "\n",
    "# View object type\n",
    "type(MeASUREs_vel)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
